
***************************************************
* Figure 4
***************************************************

use "$workdata/baseline", clear /* see $data_do/baseline.do */

replace wage2=0 if wage2==.

capture drop deflator
g deflator=.
replace deflator=wage2/7410 if year==2000
replace deflator=wage2/7711 if year==2001
replace deflator=wage2/7919 if year==2002
replace deflator=wage2/8172 if year==2003
replace deflator=wage2/8409 if year==2004
replace deflator=wage2/8577 if year==2005
replace deflator=wage2/8749 if year==2006
replace wage2=deflator

capture drop wagel*
local step=0.25
local count=-1
forvalues i=0(`step')16 {
local count=`count'+1
g wagel`count'=(wage2*(1-0.08))>`i' 
}

collapse wagel* (sum) e, by(age_q pnr)
keep if age_q>=21 & age_q<29

capture drop age_month_q_y
g age_month_q_y=floor(age_q)
tostring age_month_q_y, replace
capture drop age_month_q_m
g age_month_q_m=(round((age_q-floor(age_q))*12)+1)/100
tostring age_month_q_m, replace
capture drop age_month_q
egen age_month_q=concat(age_month_q_y age_month_q_m)
destring age_month_q, replace
g age0_q=age_q-25
g d=age0_q>=0
g dXage0_q=d*age0_q
g age0_q2=age0_q^2
g dXage0_q2=d*age0_q2
g age0_q3=age0_q^3
g dXage0_q3=d*age0_q3
g age0_q4=age0_q^4
g dXage0_q4=d*age0_q4


capture drop w e2
g w=max(0,(12/13)-abs(age0_q))
g e2=w*e

matrix x = (.,.,.,.)
preserve
local stop=16
forvalues l=0/`stop' {


local           cond="age_month_q<24.12  | age_month_q>25.04"
if `l'>=3 local cond="age_month_q<24.12  | age_month_q>25.01"
if `l'>=5 local cond="age_month_q<=25.01 | age_month_q>=25.01"
if `l'>=5 local cond="age_month_q<=25.01 | age_month_q>=25.01"
if `l'>=6 local cond="age_month_q<=25.01 | age_month_q>=25.01"

reg wagel`l' d age0_q dXage0_q  [aw=e2] if (`cond')  , vce(cluster pnr)
capture drop yhat_wage`l'
predict yhat_wage`l'
quietly su yhat_wage`l' if age0_q==0
local lo=round((r(mean)-0.05)*100)/100
local hi=round((r(mean)+0.05)*100)/100

matrix a=(`l')
matrix b=(_b[d])
matrix c=(_se[d])
matrix d=(_b[_cons])
matrix y = (a,b,c,d)
matrix x = (x\y)
}

restore


capture drop x*
svmat x
replace x1=x1/4
capture drop ci_lo ci_hi
g ci_lo=x2-1.96*x3
g ci_hi=x2+1.96*x3
capture drop nul
g nul=0 if x1!=.

twoway (rarea ci_hi ci_lo x1 if x1<=4, color(gs14)) ///
	   (line ci_lo x1 if x1<=4, lcolor(gs1) lwidth(vthin) lpattern(solid)) (line ci_hi x1 if x1<=4, lcolor(gs1) lwidth(vthin) lpattern(solid)) ///
       (connect x2 x1 if x1<=4, msymbol(o) msize(medlarge) mlcolor(gs1) mfcolor(gs8) lpattern(solid)) (line nul x1, lwidth(thin) lcolor(gs1)) ///
	   (pcarrowi -0.03 2.021272 0.01 2.021272 , mstyle(none) lcolor(gs4) lwidth(thin) lpattern(shortdash))  ///
	   (pcarrowi -0.03 2.181838 0.01 2.181838 , mstyle(none) lcolor(gs2) lwidth(vthin) lpattern(solid))  ///
	   , scheme(s1mono) legend(off) ylabel(-0.03(0.005)0.01) ytitle("RD Estimates, Employment", height(5)) ///
	   xlabel(0(0.5)4) xtitle("Earnings relative to highest SA", height(5)) name(fig4_2, replace) ///
	   subtitle("ADD Estimates over" "Earnings Distribution" "", margin(b=4) linegap(1.5))
graph save fig4_2  $figures/fig4_2, replace
graph export $figures/fig4_2.pdf, name(fig4_2) replace
graph export $figures/fig4_2.png, name(fig4_2) replace

capture drop alpha beta alpha_beta
g alpha=x4
g beta =x2
g alpha_beta=alpha+beta

twoway (rarea alpha alpha_beta x1 if alpha>=0.25 & alpha_beta>=0.25 & x1<=1.5, color(gs13))   ///
	   (rcap alpha alpha_beta x1  if alpha>=0.25 & alpha_beta>=0.25 & x1<=1.5, msize(vtiny) color(gs1))     ///
	   (line alpha x1             if alpha>=0.25                    & x1<=1.5, lcolor(gs1) lwidth(medthick) lpattern(shortdash)) ///
	   (line alpha_beta x1        if               alpha_beta>=0.25 & x1<=1.5, lcolor(gs1) lwidth(medthick) lpattern(solid)) ///
	   , legend(order(3 "< 25 years old" 4 ">= 25 years old") lcolor(gs6) lwidth(thin) col(1) position(1) ring(0) bplace(swest)) ylabel(0.35(0.05)0.55) xlabel(0(0.25)1.5) ///
	   scheme(s1mono) ytitle("Probability of earnings > x", height(5)) xtitle("Earnings relative to highest SA, x", height(5)) ///
	   name(fig4_1, replace) ///
	   subtitle("Probability of Earnings above various" "Earnings levels, just above/below age 25", margin(b=4) linegap(1.5))
graph save fig4_1  $figures/fig4_1, replace
graph export $figures/fig4_1.pdf, name(fig4_1) replace
graph export $figures/fig4_1.png, name(fig4_1) replace

graph use $figures/fig4_1
graph use $figures/fig4_2
graph combine ///
      fig4_1 fig4_2 ///
      , cols(2) rows(2) ysize(4) xsize(10) altshrink  scheme(s1mono) ///
	   name(fig4, replace)
graph save fig4  $figures/fig4, replace
graph export $figures/fig4.pdf, name(fig4) replace
graph export $figures/fig4.png, name(fig4) replace


***************************************************
* end: Figure 4
***************************************************
